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Dynamical mean field theory is employed to calculate the properties of multilayered inhomoge- 
neous devices composed of semi-infinite metallic lead layers coupled via barrier planes that are made 
from a strongly correlated material (and can be tuned through the metal-insulator Mott transition) . 
We find that the Friedel oscillations in the metallic leads are immediately frozen in and don't change 
as the thickness of the barrier increases from one to eighty planes. We also identify a generalization 
of the Thouless energy that describes the crossover from tunneling to incoherent Ohmic transport 
in the insulating barrier. We qualitatively compare the results of these self-consistent many-body 
calculations with the assumptions of non-self-consistent Landauer-based approaches to shed light 
on when such approaches are likely to yield good results for the transport. 

PACS numbers: 71.30.-|-h, 73.40.Rw, 73.20.-r, 73.40.-c 
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I. INTRODUCTION 

The fields of strongly correlated materials and of nan- 
otechnology are being united by work that investigates 
what happens when correlated materials are placed into 
inhomogeneous environments on the nanoscale. This can 
be accomplished by careful growth of strongly correlated 
materials with molecular beam epitaxy or pulsed laser 
deposition, or it may be an intrinsic property of some 
strongly correlated systems that display either nanoscale 
phase separation, or nanoscale inhomogeneity. There are 
fundamental questions about these systems — what hap- 
pens to the properties of the system when it has inhomo- 
geneities on the nanoscale and how does this spatial con- 
finement modify the quantum-mechanical correlations? 

We investigate a special case of a correlated nanostruc- 
ture, where we can carefully control the quantum confine- 
ment effects. We take a semi-infinite ballistic-metal lead 
and couple it to another semi-infinite ballistic-metal lead 
via a strongly correlated barrier material (which is from 
one to eighty atomic planes thick) . As the barrier is made 
thinner, the strongly correlated system is being confined 
in one spatial direction between the metallic leads. But 
the metallic leads induce a proximity effect on the bar- 
rier, which can deconfine the correlated system. Indeed, 
we will see that systems with a single-plane barrier still 
display upper and lower Mott bands, but they also have a 
low-energy low-weight peak to the density of states that 
arises from the proximity-effect induced states that are 
localized near the interfaces of the leads and the barrier. 
As the barrier is made thicker, this peak becomes a dip, 
which decreases exponentially with the thickness. 

We employ dynamical mean field theory (DMFT) in 
this work. This allows us to self-consistently calculate 
the properties of the inhomogeneous system, including 
Friedel-like oscillations in the leads, and the proximity- 
effect on the barrier. We do not need to make any 
assumptions about the kind of transport through this 
device, be it ballistic, diffusive, tunneling, or incoher- 
ent (via thermal excitations), since the DMFT auto- 



matically incorporates all kinds of transport within its 
formalism^. We are, however, making one approxima- 
tion in this approach — namely, we make the assumption 
that the self energy remains local, even though it can 
vary from plane to plane in the multilayered nanostruc- 
ture. Such an approximation should work fine for these 
inhomogeneous systems, since the coordination number 
remains the same throughout the device (and we are 
working in three dimensions). This is to be contrasted 
with more conventional approaches to tunneling, which 
assume a single-particle approach and employ a phe- 
nomenological potential to describe the barrier region^. 
The wavefunctions, transmission, and reflection coeffi- 
cients can be calculated, and then the transport ana- 
lyzed, as in a Landauer-based approach. In the DMFT 
calculations, we determine the potential self-consistently 
(i.e., the self energy) from the microscopic parameters of 
the Hamiltonian, and the potential can vary with the en- 
ergy of the scattering states. It is not clear that a simple 
phcnomenological potential can reproduce the same kind 
of behavior via a conventional tunneling approach. 

We assume each of the multilayer planes has transla- 
tional invariance in the perpendicular x- and y-directions. 
This allows us to use a mixed basis, Fourier transform- 
ing the two perpendicular directions to and fcj,, but 
keeping the z-direction in real space. Then for each two- 
dimensional band energy, we have a quasi one dimen- 
sional problem to solve, which has a tridiagonal represen- 
tation in real space, and can be solved with a renormal- 
ized perturbation expansion^. It is this mixed-basis rep- 
resentation (introduced by Potthoff and Nolting^) that 
allows us to solve this problem. By iterating our many- 
body equations, we can achieve a self-consistent solution. 

In addition to single-particle properties, we also eval- 
uate z-axis transport, perpendicular to the multilayers. 
Thouless introduce the idea of using the dwell time within 
the barrier to define a quantum energy scale h/tdweii, 
which turned out to describe the dynamics and transport 
of both ballistic metal and diffusive metal barriers^. 
The concept has been applied widely to the quasiclas- 
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sical theory of Josephson junctions as weli^-. If we don't 
focus on the time spent within the barrier, but instead 
try to extract an energy scale from the resistance of a de- 
vice, then we can generahze the Thouless energy to the 
case of an insulating barrier, where the transport arises 
from cither tunneling or incoherent (thermally activated) 
processes. We find that when this energy scale is on the 
order of the temperature, then we have a crossover from 
tunneling to incoherent transport. A short communica- 
tion of this work has already appeared^i. 

The organization of this paper is as follows: in Section 
II, we present a detailed derivation of the formalism and 
the numerical algorithms used to calculate properties of 
nanostructures. In Section III, we describe the single- 
particle properties, focusing on the density of states and 
the self energy. In Section IV, we generalize the concept 
of the Thouless energy, which is applied to charge trans- 
port in Section V. We end with our conclusions in Section 
VI. 



II. FORMALISM AND NUMERICAL 
ALGORITHMS 

The Hamiltonian we consider involves a hopping term 
for the electrons and an interaction term for the sites 
within the barrier region (interactions can be added in 
the metal if desired to convert the leads from a bal- 
listic metal to a diffusive metal, but we do not do 
so here). For the interaction, we employ the Falicov- 
Kimball model^ which involves an interaction between 
the conduction electrons with localized particles (thought 
of as /-electrons or charged ions) when the conduction 
electron hops onto a site occupied by the localized parti- 
cle. We consider spinless electrons here, but spin can be 
included trivially by introducing a factor of 2 into some 
of the results. The Hamiltonian is 

w = - E ^'^A^^ + E ^» U''^ - - ^) (1) 

ij i ^ / \ / 

where is a Hermitian hopping matrix, Ui is the 
Falicov-Kimball interaction, and Wi is a classical variable 
that equals one if there is a localized particle at site i and 
zero if there is no localized particle at site i (a chemical 
potential /i is employed to adjust the conduction-electron 
concentration). Since we are considering multilayered 
heterostructures, we assume that the hopping matrix is 
translationally invariant within each plane, as well as the 
Falicov-Kimball interaction. We let the z-direction de- 
note the direction where the system is allowed to have 
inhomogeneity. Then our translational invariance in the 
parameters requires that Ui — Uj if — Rj has a van- 
ishing z-component. Similarly, tij = tiiji if R^ — R^/ 
and Rj — Rj' both have a vanishing z-component, and 
Ri — Rj = Ri' — Rj' . But this requirement is quite 
modest, and allows for many complex situations to be 
considered. 



We denote the planes with a given z-component by 
a Greek label (a, /3, 7, ...). Then our requirement on 
the interaction is that Ua has a definite value for each 
plane a. The hopping matrix can have one value for 
the hopping within the plane, and different values ta,a+i 
and ta-i,Q for hopping to the plane to the right and for 
hopping to the plane to the left, respectively. For sim- 
plicity, we will only consider nearest-neighbor hopping 
here, and we assume the lattice positions R^ all lie on 
the points of a simple cubic lattice (but we do not have 
full cubic symmetry). 

Because of the translational invariance within each 
plane, we can perform a Fourier transform in the x- and 
y-coordinates to the mixed basis k^;, ky, and a (the z- 
component in real space). We define the two-dimensional 
band structure, for each plane a, by 

e2'^(k„ky) = -2ill[cosk, + cos^]. (2) 

The Green's function, in real space, is defined by 

G,,(r) = -(r.c.(r)c](0)), (3) 

for imaginary time r. The notation (O) denotes the 
trace Trexp(— — fiM])0 divided by the partition 
function Z = Trexp(— /3[7i — fiAf]) and the operators 
are expressed in the Heiscnberg representation 0{t) = 
exp(r[H - AiA/])Ocxp(-r[H - ^iN]). The symbol % 
denotes time ordering of operators, with earlier t val- 
ues appearing to the right and /3 is the inverse temper- 
ature {(3 = 1/T). We will work with the Matsubara fre- 
quency Green's functions, defined for imaginary frequen- 
cies iojn = «7rT(2n -I- 1). The Green's function at each 
Matsubara frequency is determined by a Fourier trans- 
formation 

G,,{iu;n)= / rfTe'"""G,,(r). (4) 

JQ 

We also will work with the analytic continuation of the 
time-ordered Green's functions to the real axis (retarded 
or advanced Green's functions), with iujn — > w ± iO~^. 
We use the symbol Z to denote a general variable in the 
complex plane (although we will mainly be interested in 
either Z = iw„ or Z = to + lO"*"). Finally, we work in the 
mixed basis described above, where we Fourier transform 
the X- and y-components to momentum space, to give 
GafjCk-, Z), where R^ has a z-component equal to a and 
Rj has a z-component equal to /3 (k is a two-dimensional 
wavevector). 

With all of this notation worked out, we can write the 
equation of motion for the Green's function in real spacei, 
which satisfies 

G-,\Z) = {Z + 1^)5,, - ■E,iZ)5,, + U,. (5) 

Now we go to a mixed-basis, by Fourier transforming in 
the X- and y-directions to find 

G;i(k,Z) = [Z + fi-j:^{Z)~e^''{k)]S^0 

+ taa+l5a + ll3 + taa-lSa-ll3, (6) 
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with T,a{Z) the local self energy for plane a. Finally, 
we use the identity Gaj{Z)G~^{Z) = Sap to get the 
starting point for the recursive solution to the problem: 

<5a/3 = Ga|3iKZ)[Z + ^l-^p{Z)-ef{k)] 

+ GafS-lOi, Z)tp-ii3 + Gap+liK Z)ti3+if3. (7) 



The equation of motion in Eq. |7J has a tridiagonal form 
with respect to the spatial component z, and hence it can 
be solved by employing the renormalized perturbation 
expansion'^. We illustrate the solution exactly here. The 
equation with P — a is different from the equations with 
f3 ^ a. The former is solved directly via 



G„„(k, Z) 



Go„-i(k,Z) 
G„„(k,Z) 



ta- 



(8) 



a + la 



and the latter equations can all be put into the form 

G'cta — n+l (k, Z^ta — 7i+la — n 



In a similar fashion, we define a right function and a 
recurrence relation to the right, with the right function 



G (k,Z) 

for rt > 0, with a similar result for the recurrence to the 
right. We define the left function 



La-n(^, Z) — — 



^aa— n+l(k, Z^ta—n+la—n 
Gaa-ri(k, Z) 



(10) 



and then determine the recurrence relation from Eq. ^ 



La-niKZ) = Z + ^-S„_„(Z)-ef_„(k) 

^a — na — n— it a — 71— la — n 



(11) 



We solve the recurrence relation by starting with the re- 
sult for £-oo, and then iterating Eq. up to n = 1. 
Of course we do not actually go out infinitely far in our 
calculations. We assume we have semi-infinite metallic 
leads, hence we can determine L_oo by substituting L_oo 
into both the left and right hand sides of Eq. , which 
produces a quadratic equation for L_oo that is solved by 



i-oo(k,Z) 



Z + /i-I]_oo(Z)-e2<i^(k) 



± -,/[Z + Ai-S-co(^)-^'io(k)] 



The sign in Eq. (|12() is chosen to yield an imaginary part 
less than zero for Z lying in the upper half plane, and 
vice versa for Z lying in the lower half plane. If L_oo is 
real, then we choose the root whose magnitude is larger 
than t^oo (the product of the roots equals i?.oo). In our 
calculations, we assume that the left function is equal to 
the value L-oo found in the bulk, until we are within 
thirty planes of the first interface. Then we allow those 
thirty planes to be self-consistently determined with La 
possibly changing, and we include a similar thirty planes 
on the right hand side of the last interface, terminating 
with the bulk result to the right as well. 



^9) 



(k,Z) 



(12) 



and the recurrence relation 

i?„+„(k,Z) = Z + ^-S]„+„(Z)-e^%(k) 

-Ra+n+l(ki Z) 



(14) 



We solve the right recurrence relation by starting with 
the result for Roo, and then iterating Eq. (jl4(l up to n = 
1. As before, we determine Roo by substituting i?oo into 
both the left and right hand sides of Eq. H14|l . which 
produces a quadratic equation for R^o that is solved by 



Z + A*-E^(Z)-e^(k) 



(15) 



RooiKZ) = 



± 2V[^ + /^-^-(^)-^-(k)P-4i2^. 

The sign in Eq. 1)15(1 is chosen the same way as for 
Eq. (|12|l . In our calculations, we also assume that the 
right function is equal to the value Roo found in the bulk, 
until we are within thirty planes of the first interface. 
Then we allow those thirty planes to be self-consistently 
determined with Ra possibly changing, and we include 
a similar thirty planes on the left hand side of the last 
interface, terminating with the bulk result to the left as 
well. 

Using the right and left functions, we finally obtain the 
Green's function 

GaciKZ) 



La{K Z) + RaiK Z)-[Z + fi- I]„(Z) - e, 

(16) 

where we used Eqs. ((TT|) and in Eq. ©. The local 
Green's function on each plane is then found by sum- 
ming over the two-dimensional momenta, which can be 
replaced by an integral over the two-dimensional density 
of states (DOS): 



2d 



(k)] 



GaciZ)^ J de^V'(e^')Go„(e^^Z), 



(17) 
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with 



2TrHia^ 



(18) 



and K(x) the complete elhptic mtegral of the first kind. 
If ta varies in the nanostructure, then changing variables 
to e = Ea'^/ti in Eq. ((T7|l produces 



de-^K 1 - Jl - _ G„„(4e,Z), 



(19) 

so that we can take the e variable to run from —4 to 
4 for the integration on every plane, and we just need 
to introduce the corresponding tie substitution (for e'^) 
into the left and right recurrence relations. In the bulk 
limit, where we use ta — t, we find that the local Green's 
function found from Eqs. (|17|l and (|16|l reduce to the 
well-known expressions for the three-dimensional Green's 
functions on a simple cubic lattice'^, with a hopping pa- 
rameter t. 

Once we have the local Green's function on each plane, 
we can perform the DMFT calculation to determine the 
local self energy on each plane^'^". We start with Dyson's 
equation, which defines the effective medium for each 
plane 



(20) 



The local Green's function for the ath plane satisfies 



1 



1 



GoaiZ) 



Wi 



Gll{Z)-\U' 



(21) 

with wi equal to the average filling of the localized par- 
ticles [note that this above form is slightly different from 
the usual notatiorii^, because we have made the theory 
particle-hole symmetric by the choice of the interaction 
in Eq. J^l, so that /i = corresponds to half filling in the 
barrier region and in the ballistic metal leads]. Finally, 
the self energy is found from 



Y.a{Z) = G-^l{Z)-G-\Z). 



(22) 



The full dynamical mean field theory algorithm can 
now be stated. We begin by (i) making a choice for the 
self energy on each plane. Next, we (ii) use the left and 
right recurrences in Eqs. Hll|l and (|14|l along with the 
bulk values found in Eqs. (|12|) and (|15|l and a choice for 
the number of self-consistently determined planes within 
the metal leads (which we choose to be 30 to the left 
and the right of the barrier interfaces) to calculate the 
local Green's function at each plane in the self-consistent 
region from Eqs. Hl()|l and (|19|l . Once the local Green's 
function is known for each plane, we then (iii) extract 
the effective medium for each plane from Eq. (|20|l . (iv) 
determine the new local Green's function from Eq. (|21|l . 
and (v) calculate the new self energy on each plane from 



Eq. (|22|l . Then we iterate through steps (ii)-(v) until the 
calculations have converged. 

For all of the calculations in this work, we will assume 
the hopping matrix is unchanged in the metallic leads and 
the barrier, so all ta_a±\ and all ii are equal to i, which 
we take as our energy unit. We also work at the particle- 
hole symmetric point of half filling for the conduction 
electrons and the localized electrons. This yields w\ — 
1/2 and = 0. 

There are a number of numerical details that need to 
be discussed in these computations. First, one should 
note that the recurrence relations in Eqs. Ijllfl and H14|) 
always preserve the imaginary part of i? or L during the 
recursion. Hence the recursion is stable when i? or L is 
complex. On the other hand, when they are real, we find 
that the large root is stable. Since this is the physical 
root, the recursion relations are always stable. Second, 
the integrand can have a number of singularities in it. 
When we calculate the Matsubara Green's functions, the 
only singularity comes from the logarithmic singularity 
in the two-dimensional DOS. We remove that singular- 
ity from the integration by using a midpoint rectangu- 
lar integration scheme for 0.5 < |e| < 4, and we change 
the variables for the region |e| < 0.5 from e to a;'^ = e, 
which is finite as a; — > 0, and which has a finite slope 
as a: — > 0; this allows a midpoint rectangular integration 
scheme for \x\ < (0.5)"'^/^ to accurately determine this 
second piece of the integral. When we calculate the real 
frequency Green's functions, we have the logarithmic sin- 
gularity, but we also can have a square-root singularity at 
the ath plane in the denominator of the integrand when 
ImSai^^) = and \uj + jj, ~ ReSa(a;) ~ e\ — 2. We define 
a — u) + ^ — ReI]Q(a;) + 2 and b = ui + ^ — KeJ^ai^) — 2. 
Then, if a < —4 or 6 > 4, the only singularity lies at e = 
as before. When b < —4, but —4 < a < 4, then there is 
a singularity at e = a; when a > 4, but —4 < 6 < 4, then 
there is a singularity at e = 5; and when —4 < a, 6 < 4, 
there are singularities at a and b. The singularities are 
easy to transform away by using sine and hyperbolic co- 
sine substitutions like e = u + fi — ReSa {uj) — 2 sin 9 and 
e = Lu+^—IleT,a{uj) — 2 cosh0 into the respective pieces of 
the integrands where a singularity lies. We simply deter- 
mine where all possible singularities lie (for each plane), 
set up an appropriate grid for the e variable that takes 
the different changes of integration variable into account, 
and compute the associated weight functions for the in- 
tegrations, in order to perform the integration over the 
two-dimensional DOS. Third, we find that when the cor- 
relations in the barrier are strong enough that we are in 
the Mott insulator for the bulk material, and the bar- 
rier is sufficiently thick, then the self energy develops 
a sharp structure, where the real part goes through zero 
over a small range close to w = 0, and the imaginary part 
picks up a large delta- function- like peak around co — 0. 
In order to properly pick up this behavior in the self- 
consistent solutions, we need to use a very fine integra- 
tion grid (we used up to one million points for the calcu- 
lations reported on here) to perform the integration over 
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the two-dimensional DOS. Such a fine grid is only needed 
for frequencies close to — 0, but one needs to have a 
fine enough frequency grid in lo to pick up the sharp peak 
behavior in the self energy (we use a step size of 0.001 
when there is a sharp structure in the self energy). For 
ordinary lo points, we typically used an integration grid 
of 5000 points. Fourth, these equations are easy to paral- 
lelize on the real-frequency axis, because the calculations 
for each value of frequency are completely independent 
of one another, so we simply use a master-slave approach 
and send the calculations at different frequencies to each 
of the different slaves until all frequencies are calculated. 
This approach has an almost linear scale up in the par- 
allelization speed. 

In addition to these single-particle properties, we also 
are interested in transport along the z-axis (perpendic- 
ular to the multilayered planes). The resistance of the 
nanostructures can be calculated by a Kubo-based linear 
response formalismpii (i.e., a current-current correlation 
function). We begin with the current operator at the ath 
plane 



Jz — ^ ' jzct; 



i in 2d plane 

This operator sums all of the current flowing from the 
ath plane to the a -\- 1st plane. 

The current-current correlation function is defined to 

be 

n^piii^i)^ r dTe^'"^%3Ur)iz0m, (24) 



with ii^i = iTrT2l the Bosonic Matsubara frequency and 
with the dc conductivity matrix determined by the ana- 
lytic continuation of Eq. (|24|l to the real frequency axis 
via 



cTapiv) = lim Re 



(25) 



Substituting Eq. (|23|l into Eq. 124|) . evaluating the con- 
tractions in terms of the single-particle Green's func- 
tions, performing the integration over r to convert to the 
Matsubara frequency representation, and performing a 
Fourier transform over the 2d-spatial coordinates, yields 
the following result after some straightforward algebra: 



G/3+ia+i(k, zcj„i)G'a/3(k, iuj„i + ivi) 
G/3Q+i(k, iw,„)Ga/3+i(k, iuj^ + m) 



Now we need to perform the analytic continuation from 
the imaginary to the real frequency axis^'^. This is done 
by first converting the summations over the Matsubara 
frequencies into contour integrals that enclose all of the 
Matsubara frequencies and are multiplied by the Fcrmi- 
Dirac distribution function f{uj) = l/[l + cxp(/3cj)] which 
has a pole at each Matsubara frequency. Then the con- 
tours are deformed to go along lines parallel (but just 
above or just below) the real axis, and the real axis 
shifted by —ivi- At this point we replace f{uj — ivi) by 
/((jj) and then analytically continue ivi v + zO^. The 
algebra is once again straightforward but somewhat te- 
dious. The final result is 



k L 

Gai3 (k, uj + iy)lniG(3+ia+i (k, uj) 
+ G'Q/3+i(k,a; -I- zy)ImG/3a+i(k,cj) 
-I- GQ+i/3(k,a; + z^)ImG,3+iQ(k,cj) 

- GQ+i/3+i(k, w + i^)ImG;3a(k, w)| 

- G'^+iQ+i(k,t^)ImGa/3(k, w + i^) 
Gfl„+i(k, w)ImGa/3+i(k, w + v) 



+ 
+ 



GWi^ik, w)ImGa+i/3(k, uj - 



G^„(k,a;)ImGa+i/3+i(k,Lj + 1/)| 



(27) 



The last step is to evaluate the dc conductivity matrix, 
which becomes 

ImG/3a+i (e, Lu)lu\Gai3+i (e, w) 
-I- ImG/3+ia(e, w)ImGQ+i/3(e, w) 
- ImG^+ia+i(e,Lj)ImGa/3(e,^) 
+ ImG/3a(e,Lj)ImGa+i/3+i(e,cj) . (28) 

The conductivity matrix has the dimensions e^/ha^, 
which is the inverse of the resistance unit, divided by 
two factors of length, and is the correct units for the 
conductivity matrix. 

Since the conductivity matrix is not as familiar as the 
scalar conductivity used for homogeneous problems, we 
will briefly derive how one extracts the resistance of the 
nanostructure from the conductivity matrix. The key 
element that we use is that the current density that flows 
through each plane is conserved, because charge current 
can neither be created nor destroyed in our device. The 
continuity equation, then says that the current density 
through the ath plane, /q,, is related to the electric field, 
Ej3, between the /?th and (3 -\- 1st plane via 



Ia= crai3{0)Ei3 = /, 



(29) 
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where we set the current density on each plane equal to 
a constant value I. Inverting this relation to determine 
the electric field gives 



(30) 



The voltage across the nanostructure is just the sum of 
the electric field between each plane, multiplied by the 
interplane distance (we assume a constant dielectric con- 
stant throughout), so we can immediately determine the 
resistance-area product (specific resistance) from Ohm's 
law 



V 



RnCL = — = } ya 

a/3 



(31) 



One needs to pursue a similar type of analysis to examine 
the thermal transport properties (thermopower and ther- 
mal resistance), but it is somewhat more complicated, 
because the thermal current is not conserved from one 
plane to another plane, as is the charge current. We will 
present results for such a calculation elsewhere (at half 
filling, where we restrict ourselves in this paper, there is 
no thermopower by particle- hole symmetry). 

The only mathematical issue associated with this anal- 
ysis is that we have assumed the conductivity matrix is 
invertible. In general, this is not true when there is no 
scattering in the metallic leads. In this case, we need 
to truncate the conductivity matrix to consider only the 
block that covers all of the planes in the barrier and the 
first metallic plane to the left and to the right of the 
barrier. This matrix is always invertible, and allows cal- 
culations to be performed easily (if we were to include 
a larger matrix, we find that the resistance does not in- 
crease as we increase the number of planes within the 
metallic leads that we include in the conductivity matrix 
block that is inverted, at least until we run into precision 
issues for the calculations). Of course, if the metallic 
leads have scattering, there are no numerical issues asso- 
ciated with the matrix inversion (except when the matrix 
is made too large and the system has approached the bulk 
limit, see below), but we need to decide how far down the 
metallic leads we will perform the actual measurement, 
since the voltage grows with the thickness of the metallic 
leads included in the calculation (when there is scattering 
in the leads). 

In order to calculate the dc conductivity matrix in 
Eq. (|28|) . we need to evaluate the off-diagonal compo- 
nents of the Green's functions. This is easy to do using 
the renormalized perturbation expansion, and the right 
and left functions. We find two recurrence relations 



Gaa-n{^, ^) — 

(defined for n > 0) and 



GctOL-\-n 



G'aa — n+l^a — n+la — n 

La-„{e,uj) 



^aa+n — l^a+n— la+n 



(32) 



(33) 



(also defined for n > 0). The other off-diagonal 
Green's functions are found from the symmetry relations: 

Gcta — n — Ga — na and Ccta+n — Ga-^-nQ,- 

The computation of the junction resistance for a given 
temperature is relatively simple to perform. First, one 
must calculate all of the local self energies for each plane, 
using the algorithm described above. Then, for each fre- 
quency Lo, one can calculate all of the Green's functions 
that enter into the formula for 0-0/3(0). It is best to eval- 
uate the integral over uj for many different temperatures 
"at the same time" since the only thing that changes with 
temperature (when at half filling, where the chemical po- 
tential is fixed and does not vary with T) is the Fermi 
factor derivative. Since evaluating at each frequency is 
independent of every other frequency, this algorithm is 
also "embarrassingly parallel" . 

One final comment is in order about the formalism for 
calculating the junction resistance. Namely, how does 
it relate to a Landauer approach to the resistance? In 
the Landauer approach^ one does not calculate a con- 
ductivity matrix, but instead determines the transport 
directly by evaluating the Green's function Gap where a 
lies at the left interface and /? lies at the right interface. 
We believe one can show that these two approaches are 
completely equivalent if one uses the same self energies 
for the inhomogeneous structure to calculate the Green's 
functions that enter into the transport calculation. We 
will examine this relationship in a future publication. 

In a homogeneous (bulk) noninteracting system, we 
find that the Green's functions satisfy 



Gaa^n 



^4^2 ~{u; + fi- e)2 
ti) + /i — e 



(34) 



+ 1 



. v/4t2 - (co + /i - e)2 



when e lies within the band [|a; + /i — e| < 2]. Note 
that ImGQ./3(e, cj) is not always negative when a ^ (3. 
This occurs because we are using a mixed basis, and the 
imaginary part of the Green's function does not have 
a definite sign in this basis. We can substitute these 
Green's functions into the expression for the conductivity 
matrix, to evaluate the result for the bulk. We find the 
matrix has all of its matrix elements equal to each other, 
and they assume the value 



O'a/3(0) 



e 



dep^''{e) « 0.63 



/ia2' 



(35) 



for the case of half filling fi = (since every matrix 
element is the same, the conductivity matrix is not 
invertible, but the resistance can still be calculated). 
This result will lead to precisely the Sharvin contact 
resistancoi^ii^ii^ when we convert the conductivity into 
a resistance (the resistivity of a ballistic metal vanishes, 
but the resistance is nonzero). 
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III. SINGLE-PARTICLE PROPERTIES 

We perform our calculations at half filling {jj, ~ 0, 
(cj^Cj) = 1/2, and wi = {wi) = 1/2). This has a num- 
ber of advantages. First, because the chemical poten- 
tial is the same for the metallic leads and the barrier, 
there is no electrochemical force that reorganizes the 
electrons to a screened dipole layer at each of the inter- 
faces, instead the filling remains homogeneous through- 
out the system. Second, the chemical potential is fixed 
as a function of temperature, so there is no need to per- 
form imaginary-axis calculations to determine the chem- 
ical potential as a function of temperature. We usually 
calculate the Matsubara Green's functions anyway, to 
test the accuracy of the real-axis Green's function, by 
comparing the Matsubara Green's functions calculated 
directly with those calculated from the spectral formula 
via the real-axis DOS (usually the accuracy is better than 
three decimal points for every Matsubara frequency). 
Third, we can perform calculations of the resistance at all 
temperatures in parallel, because the chemical potential 
does not vary with temperature (recall, the DOS of the 
Falicov-Kimball model is temperature independent for 
the DMFT solutior.^^). Fourth, the particle-hole symme- 
try of the DOS allows us to have another check on the ac- 
curacy of the calculations because we do not invoke that 
symmetry in our calculations. Fifth, there is a metal- 
insulator transition (MIT) in the bulk Falicov-Kimball 
model on a cubic lattice when U « 4.9i, so the solutions 
at half filling include the MIT. For these reasons, we find 
this case to be the simplest one to consider in a first 
approach to the inhomogeneous many-body problem. 

We also reduce the number of parameters in our cal- 
culations by assuming all of the hopping matrix elements 
are equal to t for nearest neighbors. This is by no means 
necessary, but it allows us to reduce the number of pa- 
rameters that we vary in our calculations, which allows 
us to focus on the physical properties with fewer calcu- 
lations. The hopping scale t is used as our energy scale. 
We also include 30 self-consistent planes in the metallic 
leads to the left and to the right of our barrier, which is 
varied between 1 and 80 planes in our calculations. 

The first problem we investigate is the extreme quan- 
tum limit of having one atomic plane in the barrier of our 
device. We tune the Falicov-Kimball interaction in the 
one barrier plane from U — 1 to U — 20, which goes from 
a dirty metal to well into the Mott insulating regime. But 
the Mott insulating phase does not like being confined to 
a single atomic plane, and there is a metallic proximity 
effect, where the metallic DOS leaks into the insulator 
DOS at low energies. The result is that we do not expect 
the single-plane barrier to be too resistive. This is easiest 
to see when we consider the local DOS within the barrier 
plane, as plotted in Fig. ^ There we see that the DOS 
starts to be reduced at the chemical potential as we in 
crease U, but there is still substantial DOS at the Fermi 
energy when U ~ 4.9. In fact, as U is increased, we 
see that the upper and lower Mott-Hubbard bands form. 
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FIG. 1: Barrier DOS as a function of the Falicov-Kimball 
interaction U. The different line widths and styles denote 
different U values, as detailed in the legend. Note how the 
DOS initially evolves as in the bulk, with the DOS being 
reduced near a; = 0, and the band width increasing. But as 
we pass through the Mott transition, we see that the double- 
peak Mott-Hubbard bands appear, but so does a low-energy 
(interface-localized) band near uj = 0, which looks like a low- 
weight metallic band for large U. 



centered at zLU/2, but there is significant DOS that re- 
mains centered at a; = 0, and it even develops a small 
peak for U > 10. The origin of, and the size of this 
peak, can be shown to arise naturally from the renor- 
malized perturbation theory expressions for the Green's 
functions, but we do not do so herein. We anticipate that 
these states are localized at the interface, and represent 
the states that an incident electron can tunnel through 
to go from one metallic lead to the other in a transport 
experiment. These results show a number of interesting 
features of the coupling of a Mott insulator to a metallic 
lead: (i) the Mott transition remains in the sense that 
Mott-Hubbard bands continue to form, with their ori- 
gin clearly seen near the MIT; (ii) the interface-localized 
states have a metallic character (i.e., a peak at a; = 0) in 
the large-C/ regime; and (iii) the proximity effect appears 
to always be active, and able to create states within the 
barrier at low energy, but the total weight in those states 
is low, so medium to high energy properties of the Mott 
insulator phase will remain similar to the bulk. 

Next we examine what happens as we increase the 
barrier thickness for given values of U. Our focus is 
on three generic values of interest: U = 2, which is a 
strongly scattering, diffusive metal; J7 = 4, which is so 
close to the MIT, that the bulk DOS show a significant 
dip near uj — 0; and U = 6, which is well within the 
Mott-insulating phase. We first examine how the metal- 
lic leads are influenced by the presence of the barrier. 
We set the origin of the a variables so that a = cor- 
responds to the first barrier plane (hence planes —1 to 
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FIG. 2: Lead DOS for an TV = 5 barrier device with U = 2. 
The different panels show the DOS in the first metal plane to 
the left of the barrier, in the second, the third, the tenth and 
the thirtieth. Note how the system approaches the bulk cubic 
DOS as it moves further from the interface, as expected. A 
careful examination of the panels shows that the "flat" region 
with \lo\ < 2 shows a half-period oscillation for each unit 
of distance from the current plane to the interface, but the 
amplitude shrinks dramatically as we move further from the 
interface. 



—30 represent the thirty planes to the left of the bar- 
rier, with —1 closest to the barrier). In Fig. 12 we show 
results for U = 2 and five representative planes in the 
metal (the device has five barrier planes) . In Fig. 13 we 
show the same results for L'^ = 4 and in Fig. 0] we show 
the same results for U = 6. The first thing to notice is 
that the DOS is close to that of the bulk simple cubic 
lattice for 30 planes away from the interface, indicating 
that our choice of thirty self-consistent planes is reason- 
able. Next, note that the amplitude of the oscillations 
grows as U increases. Third, the number of half peri- 
ods in the oscillations increases with the distance away 
from the interface (both for < 2 and |a;| > 2). The 
source of these oscillations is the Friedel oscillations (with 
a wavelength on the order of two lattice spacings for half 
filling) that we expect associated with the disturbance 
of the Fermi sea of the metal by the proximity to the 
interface. 

There are two interesting questions to ask about these 
results: how thick does the barrier have to be before the 
Friedel oscillations become frozen in the metallic leads 



FIG. 3: Lead DOS for an TV = 5 barrier device with [/ = 4. 
The different panels show the DOS in the flrst metal plane 
to the left of the barrier, in the second, the third, the tenth 
and the thirtieth. Note how the amplitude of the oscillations 
increases as U increases. 



and don't change with a thicker barrier, and do we see 
oscillatory behavior in the barrier, where we instead ex- 
pect there to be exponentially decaying wavefunctions? 
We find that the answer to the first question is that the 
structure is already essentially frozen in for a single-plane 
barrier, and it does not evolve much with the barrier 
thickness (although it does show much evolution with 
the interaction strength). This perhaps sheds some light 
on why non-self-consistent Landauer based approaches 
for transport have been so successful. If one has a good 
guess for the semi-infinite lead DOS, then it does not 
change much as the thickness increases, so that guess 
will work well for all calculations with the same strength 
of electron correlations. 

To examine the second question, we plot results for the 
DOS at a fixed frequency (four chosen for each U value) 
in Fig. There are six different thicknesses plotted for 
each U value. The curves all lie on top of each other 
for the metallic lead planes, indicating that the Friedel 
oscillation structure is frozen in starting at TV = 1 (and 
we can read off the oscillation wavelength to be two lat- 
tice spacings, with a sharp decrease of the amplitude as 
one moves away from the interface). In the barrier, we 
see that there are only oscillations close to the interface, 
then the curves either flatten out or exponentially decay 
with thickness. But the curves continue to lie on top of 
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FIG. 4: Lead DOS for an TV = 5 barrier device with [/ = 6. 
The different panels show the DOS in the first metal plane to 
the left of the barrier, in the second, the third, the tenth and 
the thirtieth. Note how the amplitude of the oscillations is 
even larger here. A careful examination shows there are also 
oscillations (with the same kind of increase in the number of 
half periods with the distance from the interface) in the region 
\uj\ > 2. 



each other (except for the middle plane of the barrier 
for small uj and U = 6). These results, once again, show 
that another of the assumptions of the non-self-consistent 
Landauer-based approaches, that there is an exponential 
decay with a well defined decay length in the insulat- 
ing barrier regions, holds here as well, but one needs to 
properly predict the decay length to perform accurate 
calculations. 

Our final summary of the DOS is included in false color 
plots (the color, or grayscale, denoting the height of the 
DOS at a given plane) to emphasize the spatial location 
and amplitudes in the oscillations. Fig. |B1 shows the re- 
sults for = 1 and U ~ 6 and Fig. [7| shows the results 
with = 20 and U — 6 (only half of the nanostruc- 
ture planes are shown due to the mirror symmetry) . The 
color scale (or grayscale) needs to use a banded rainbow, 
with the different colors (grayscales) separated by bands 
of black in order to pick up the small amplitude oscilla- 
tions in the background of the large DOS. Note how the 
Friedel oscillations are essentially the same in the two 
plots, indicating this freezing of the oscillations starting 
at = 1. There are also oscillations visible near the 
metal band edges, indicating Friedel-like oscillations due 
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FIG. 5: DOS at specific values of a; as a function of the plane 
position in the device. We plot only the left-hand piece of the 
plots, since the right-hand piece is a mirror image of the left- 
hand piece. Note that the f/ = 6 panel is a semilogarithmic 
plot. The four values of for U = 2 are 0.0, 3.0, 4.0, and 
5.0. The barrier thicknesses are = 1, 5, 10, 20, 40, and 
80. The four values of tj for (7 = 4 are 0.0, 2.5, 3.5, and 5.0. 
The barrier thicknesses are A' = 1, 5, 10, 20, 40, and 80. The 
four values of for [/ = 6 are 0.0, 0.2, 0.4, and 1.0. The 
thicknesses are A' = 1, 4, 7, 10, 15, and 20. Note how all 
curves lie on top of each other in the metallic lead, indicating 
the structure in the metallic lead is frozen in for an A'^ = 1 
barrier, and does not significantly change with increasing N . 
In the barrier, we only have oscillations at the interface, and 
then the curves either are flat with thickness il] — 2 and 4), 
or exponentially decreasing or flat {U = 6). The little tails 
that stick out for the lowest two frequencies with 17 = 6 show 
that the middle plane of the barrier does not follow the same 
exponential decay as the other planes do. But the exponent 
of the exponential decay is frozen in starting at A ~ 1. 



to the different total bandwidths of the two materials 
joined in the nanostructure. The DOS in the barrier at 
low frequency becomes very small very quickly on these 
linear scales, but it is nonzero (see Fig. EJ. 

The final single-particle property we consider is the 
imaginary part of the self energy at the central plane of 
the barrier at low energy in Fig.lHI In the bulk, the imag- 
inary part of the self energy vanishes within the Mott- 
Hubbard gap, except for a delta function at u; = whose 
weight can be used as a quasi-order parameter for the 
Mott transition at half filling (but not away from half 
filling—). In the nanostructures, the imaginary part of 
the self energy never vanishes in the bulk gap region, but 
it can assume very small values, with a sharp peak, of 
finite width, developing at = 0. This peak grows in 
height and narrows as the barrier is made thicker. It is 
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FIG. 6: False-color plot of the DOS for a TV = 1 barrier plane 
device with U = 6. The barrier plane is just the lowest plane 
at the bottom of the figure, while the thirty metallic planes 
lie on top. Note how the ripples of the Friedel oscillations 
are most visible in the central region, where the DOS has a 
plateau. (Color version online.) 
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FIG. 7: False-color plot of the DOS for a iV = 20 barrier 
plane device with U — 6. The barrier planes are the lower 
ten planes, while the thirty metallic planes lie on top. Note 
how the ripples of the Friedel oscillations agree with those 
in Fig. |S| In the barrier, the DOS decreases rapidly on this 
linear scale, and shows few oscillations, but one can see some 
small oscillations near the band edges in both regions. (Color 
version online.) 
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FIG. 8: Semilogarithmic plot of the imaginary part of the self 
energy on the central plane of the barrier at small frequency 
for five different thickness barriers (A'^ = 1, 4, 7, 10, and 
15). Note how the imaginary part of the self energy becomes 
very small for frequencies close to oj = 0, but as we approach 
u — 0, sl sharp delta- function-like peak develops that narrows 
as the barrier is made thicker. It is precisely this structure 
that is hard to reproduce with numerical calculations. Note 
that this kind of a self energy is very similar to what is seen 
in the hypercubic lattice in infinite dimensions. 



point whether such behavior could lead to enhancements 
in the nanostructures, even though the self energy has 
similar properties. 



IV. GENERALIZED THOULESS ENERGY 



a challenge to try to calculate such a structure numeri- 
cally, especially due to the loss of precision in extracting 
the self energy from the Dyson equation during the itera- 
tive algorithm. It requires a fine enough frequency grid to 
pick up the narrow structure, and it requires a sufficiently 
fine integration grid for e, in order to accurately deter- 
mine the peak value. Note how the self energy evolves 
from a relatively broad featureless structure to a very 
sharply peaked structure as the barrier is made thicker. 
This kind of a peaked self energy is similar to what is 
seen in the exact solution on the hypercubic lattice in in- 
finite dimensions. There the Mott transition is actually 
to a pseudogap phase, with the DOS vanishing only at 
the chemical potential, but there is a region of exponen- 
tially small DOS in the "gap region" . The sharp features 
in the self energy led to a significant enhancement of the 
low-temperature thermopower on the hypercubic lattice, 
when the system was doped off of half fiUingiS (and wi 
changed to produce an insulator). It is unclear at this 



It is important to try to bring semiclassical ideas 
of transport into transport in nanostructures, to see 
whether those concepts have useful quantum analogues. 
Thouless was the first to investigate such ideas for diffu- 
sive metal barriers^. He considered the idea of a dwell 
time in the barrier for an electron that tries to travel 
through the barrier. If we assume the electron takes a 
random walk through the barrier, then the time it spends 
inside the barrier is proportional to the square of the 
thickness of the barrier (with the proportionality being 
related to the diffusion constant). Since one can extract 
the diffusion constant, via an Einstein relation, from the 
junction resistance, Thouless could construct a quantum- 
mechanical energy h/tdweii from these classical ideas. It 
turns out that this energy scale plays a significant role 
in determining the quantum dynamics of many different 
kinds of nanostructures. For example, it can be easily 
generalized to take into account ballistic metals, where 
tdweii = Na/vp for a barrier of thickness Na, with vp 
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the Fermi velocity. The Thouless energy appears to be 
the critical quantum energy scale that determines the 
dynamics through weakly correlated nanostructures; its 
success in the theory of Josephson junctions is particu- 
larly noteworthy^. 
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FIG. 9: Thouless energy for a. U — 4 (diffusive, but very 
strongly scattering metal) barrier as a function of the barrier 
thickness L = Na. The different curves correspond to dif- 
ferent temperatures. The top panel multiplies the Thouless 
energy by to try to isolate the prefactor for the diffusive 
transport, while the bottom panel plots the Thouless energy 
on a semi-logarithmic plot. Note that the temperature de- 
pendence of the constant, seen for thick barriers in panel (a), 
arises from the fact that the U = 4 DOS has significant low- 
energy structure, because there is a dip that develops near the 
chemical potential, so the temperature dependence is both 
stronger than expected for normal metals, and anomalous be- 
cause many more states are involved as T is increased, i.e. it 
behaves more like an insulator. 

So the fundamental question we wish to investigate is 
can the concept of a Thouless energy be generalized to 
a strongly correlated system, where transport through 
a nanostructure is either via tunneling or via incoher- 
ent thermal excitation. The answer is yes, and we do so 
by first trying to extract an energy scale from the resis- 
tance of the junction, which is able to track the putative 
thermal dependence of the resistance when we are in the 
incoherent thermal transport regime. A simple dimen- 
sionality argument shows that the form 



Th 



i?„a22e2j duj[~df /duj]pbuik{^)Na 



(36) 



has the the kind of dependence we are looking for. The 
symbol Pbuik{'-^) is the local DOS in the bulk for the ma- 
terial that sits in the barrier of the nanostructure. If 
we check the dimensions, we see that i?„ has dimensions 
h/e^, and the DOS has dimensions so Eth is an 



energy [note Eq. (|36|) corrects typos in an earlier worki] . 
When we examine systems where the barrier is a metal, 
then at low temperature the bulk DOS can be replaced by 
a constant in the integral, and we reproduce the known 
forms for the Thouless energy for ballistic [Eth ~ C /Na) 
and diffusive {Eth ~ C' /[NaY) electrons because the re- 
sistance is independent of the thickness for a ballistic 
metal barrier and it grows linearly with the thickness for 
a diffusive metal barrier. This method of generalizing the 
Thouless energy also avoids us having to try to answer 
the question of how long does it take an electron to tun- 
nel from the left to the right lead, and it reproduces all 
of the known forms for the Thouless energy in a unifying 
formula that does not require us to even use the Einstein 
relation to extract a diffusion constant or to determine 
the Fermi velocity for an anisotropic Fermi surface (in 
the ballistic case). 

We plot the results for this Thouless energy as a func- 
tion of thickness in Fig. for [/ = 4. In panel (a), we 
multiply Eth by the square of the length L = Na of 
the barrier. The different curves correspond to differ- 
ent temperatures. If the Thouless energy went exactly 
like C'/L^, then all of the curves would be straight lines, 
with a temperature-dependent value C'(T). But we see 
some curvature for small barrier thicknesses. This arises 
mainly from the fact that in addition to the diffusive con- 
tribution to the resistance, there is a contact resistance, 
so for thin barriers, we do not have a pure 1/L^ behavior. 
Note, however, that the Thouless energy has little tem- 
perature dependence at low temperature, as expected. In 
panel (b), we plot the curves on a semi- logarithmic plot, 
so one can see how small the Thouless energy becomes 
for thicker junctions. 

The Thouless energy is plotted versus temperature on 
a log- log plot for J7 = 6, which corresponds to a Mott- 
insulating barrier with a small correlation-induced gap. 
The dashed line indicates where Eth — T, which is an 
important crossover point for dynamics, as we will see 
below. Note that the temperature dependence is sig- 
nificant in an insulator, because the integral in the de- 
nominator of Eq. I|3t)|) has strong temperature depen- 
dence in the insulator, but the resistance does not in 
the tunneling regime at low temperature. If we used 
the Thouless energy to determine the tunneling time via 
ttunnei — fi/ Exh, wc would find tunneling times rapidly 
approaching zero asT —^ 0. We will not comment further 
here as to whether there is any substance to using such 
results to describe the quantum dynamics of the tunnel- 
ing process. Instead we simply want to conclude that 
the concept of the Thouless energy can be generalized to 
strongly correlated systems, and we will see below that 
the crossover point where Eth ~ T has important physi- 
cal interpretations that will be developed in the next sec- 
tion. Finally, the generalization of the Thouless energy 
to correlated systems changes the idea of a single energy 
scale being associated with the transport, since now the 
energy scale develops strong temperature dependence. If 
a single number is desired, then we would propose to use 
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FIG. 10: Thouless energy for a U — 6 (Mott-insulating) 
barrier as a function of temperature on a log-log plot. The 
different curves correspond to different thicknesses of the bar- 
rier, ranging from A'' = 1 for the top curve to A'^ = 2, 3, 4, 
5, 7, 10, 15, and 20 as we move down the plot. Note how 
the Thouless energy picks up dramatic temperature depen- 
dence here. The dashed line is the curve where Eth, = T. We 
find that when the Thouless energy equals the temperature, 
interesting effects occur (see below). 



the energy scale where the Thouless energy is equal to 
the temperature, indicated by the points of intersection 
of the solid lines with the dashed curves in Fig. ^| 



V. CHARGE TRANSPORT 

The dc resistance is a low-energy property of the nanos- 
tructure, and so it requires the results of the single- 
particle properties to be determined accurately at low 
energy. This is not difhcult for metallic barriers with any 
degree of scattering, as long as the numerical subtleties 
discussed above are taken into account in the analysis, 
but it does create problems for thick Mott insulators. 
We need to be able to properly determine the structure 
seen in Fig.|51as the barrier is made thicker, and this can 
exhaust the numerical resources, or the numerical preci- 
sion available for a given calculation. For our work, we 
were not successful in examining U — 6 barriers thicker 
than N ^20. 

We plot the resistance-area product in Fig. ^] for 
T = 0.01 and four different U values: C/ = 2, a diffusive 
metal near the loffe-Regel limit of a mean free path on 
the order of a lattice spacing; U — 4, a strongly scatter- 
ing, anomalous metal, that has a strong dip in the DOS 
near the chemical potential; [/ = 5, a Mott-insulator 
that is nearly critical; and U — 6, a, Mott-insulator with 
a small correlation- induced gap. In panel (a), we have 
a semi- logarithmic plot, which is useful for picking out 
tunneling behavior via an exponential increase of the re- 
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FIG. 11: Resistance-area product for nanostructures with 
U = 2, 4, 5, and 6, and various thicknesses. Panel (a) is a 
semi-logarithmic plot, while panel (b) is a linear plot. The 
temperature is T = 0.01 in both panels. Note how the cor- 
related insulator {U = 6) has an exponential growth with 
thickness as expected for a tunneling process, but it turns 
over at the thickest junction, indicating a crossover to the in- 
coherent transport regime. The U — 5 data, which is close to 
the critical point for a MIT, has neither linear, nor exponen- 
tial growth of its resistance-area product. The metallic cases 
{U = 2 and 4) have perfect linear scaling of the resistance 
with current, with a nonzero intercept corresponding to the 
contact resistance. This may be surprising for (7 = 4, be- 
cause it is so strongly scattering (with a mean free path much 
less than a lattice spacing), that one would not think a semi- 
classical approach should apply there. The constant satisfies 
CTo = 2e^/ha'^. 



sistance with thickness. This is clearly seen for the Mott 
insulator with J7 = 6, with the beginnings of a crossover 
occurring near N — 20, but the near-critical insulator 
at U = 5 does not grow exponentially, nor does it grow 
linearly [see panel (b)]. The data for U — 2 and U — 4, 
both show linear increases with thickness, with a nonzero 
intercept on the y-axis denoting the nonzero contact re- 
sistance with the metallic leads. It is surprising that 
this linear "Ohmic" scaling holds for systems that are so 
strongly scattering, that their mean free path is much 
less than one lattice spacing. 

Our final figure plots the resistance-area versus tem- 
perature for {a) U = 4 and (b) f/ = 6 [Fig. In 
panel (a), we can infer a linear dependence of -R„a^ ver- 
sus L for all temperatures, so this barrier is always Ohmic 
in nature. But it has quite anomalous temperature de- 
pendence, looking like an insulator, whose resistance is 
reduced as the temperature increases. In panel (b), we 
see an exponential dependence of RnO,^ versus L at low 
temperature, marked by the equidistant step increases of 
RnO,^ as the thickness increases (recall this is a log- log 
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This is the incoherent "Ohmic" regime for the transport. 
The sohd dots represent the resistance-area product at 
the Thouless energy, determined by finding the temper- 
ature where Eth = T from Fig. 1101 and marking those 
points on the curves in panel (b) . A dashed hne guide to 
the eye is drawn through these points. One can clearly 
see that the point where the Thouless energy equals the 
temperature determines the crossover from tunneling to 
incoherent transport. Surprisingly, this crossover occurs 
at a lower temperature for a thicker barrier. This occurs, 
because the tunneling resistance is higher for a thicker 
barrier. As T increases, the Ohmic resistance, deter- 
mined by multiplying the temperature-dependent bulk 
resistivity by the thickness and dividing by the area, will 
decrease. Once it is essentially equal to the tunneling re- 
sistance, there will be a crossover from tunneling, which 
provides a "quantum short" across the junction for low 
T, to "Ohmic" (incoherent) thermally activated trans- 
port. This must occur at a lower temperature for more 
resistive junctions, and hence the thicker junctions have 
the crossover before the thinner junctions. Note that the 
temperature scale for this crossover does not appear to 
have any simple relation to the energy gap of the bulk 
material, instead it is intimately related to the dynami- 
cal information encoded in the generalized Exh found in 
Eq. 123). 

Wc do not consider thermal transport there, since the 
thcrmopower vanishes for this particle-hole symmetric 
case and the thermal resistance is not as interesting in 
systems with vanishing thermopower. 



FIG. 12: Resistance-area product for nanostructures with (a) 
U = A and (b) ?7 = 6 as a function of temperature [panel (a) 
is on a linear scale, and panel (b) is a log-log plot]. In panel 
(a) we include results for A'^ = 1, 2, (lowest two curves), 5, 10, 
15, 20, 40, 60, and 80. Note how at each temperature there 
is a linear dependence of the resistance-area product with the 
thickness of the junction. Note further, that these junctions 
have anomalous temperature dependence for a metal (they ac- 
tually look insulating in their dependence). In panel (b), we 
show the results for U = Q with AT = 1 - 10, 15 and 20. Note 
at low temperature we have tunneling, as the resistance-area 
product is weakly dependent on temperature, and the steps 
are equally spaced as a function of thickness, indicating expo- 
nential dependence on the thickness. At higher temperatures, 
there is a crossover to the incoherent transport regime, with 
the resistance-area product picking up a strong T dependence, 
and scaling linearly with the thickness. The dotted line that 
connects the solid dots is a plot of the resistance-area value 
at the temperature where Exh = T which determines the 
crossover. 



plot). The temperature dependence is also weak in this 
region, indicated by the flatness of the curves. Hence the 
system is in the tunneling regime at low temperature. 
As T rises, there is a relatively sharp crossover region, 
where i?„a^ begins to pick up strong (exponentially ac- 
tivated) T dependence, and i?„a^ grows linearly with L. 



VI. CONCLUSIONS 

In this contribution we worked with a generalization of 
DMFT to inhomogeneous systems to calculate the self- 
consistent many-body solutions for multilayered nanos- 
tructures that have barriers that can be tuned to go 
through the Mott transition. We developed the compu- 
tational formalism thoroughly (based on the algorithm of 
Potthoff and Nolting), and although we applied it only 
to the Falicov-Kimball model, it is obvious that one can 
trivially add mean-field-like interactions such as Zeeman 
splitting for magnetic systems, or long-range Coulomb 
interactions for systems with mismatched chemical po- 
tentials. In addition, one can invoke whatever impurity 
solver desired for the local DMFT problem on each plane, 
which extracts a new self energy from the current local 
Green's function. We studied both the single-particle 
properties and the charge transport. 

There are a number of interesting results that came 
out of this analysis. First, we found that as the strength 
of the correlations increases in the barrier, there is a 
stronger feedback effect on the Friedel-like oscillations 
that appear in the metallic leads, but those oscillations 
vary little with the thickness of the barrier for a fixed 
interaction strength. Second, there are few oscillations 
inside the barrier except close to the interface with the 
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metallic leads, but the behavior in the barrier, of either 
an exponential decay, or of a constant DOS, gets frozen 
in for a relatively thin barrier, and the DOS changes lit- 
tle with increasing the thickness of the barrier, except 
when there is exponential decay which will always de- 
crease within the correlation-induced gap. Third, the 
Mott insulating barrier develops a narrow peak-like struc- 
ture in the imaginary part of the self energy that ap- 
proaches the bulk delta function result. This narrow and 
tall peak is difficult to determine accurately with the nu- 
merics and limits the ability to study thick insulating 
barriers. Fourth, we showed how to generalize the con- 
cept of a Thouless energy to become a function of T for a 
strongly correlated Mott insulator. Our unifying form for 
the Thouless energy includes the results for both the bal- 
listic and diffusive metals as well. We identified an energy 
scale that describes the crossover from tunneling to inco- 
herent transport in these nanostructures; it corresponds 
to Exh = T. This energy scale is quite useful in other 
areas such as in the theory of Josephson junctions, which 
will be presented elsewhere. Sixth, we analyzed the resis- 
tance of these devices and found interesting behavior, in- 
cluding anomalous metallic behavior (but no tunneling) 
for a strongly scattering metal, and the crossover from 
tunneling to Ohmic transport for insulating barriers. 

This work also shed light on other approaches to trans- 
port through multilayered structures like the Landauer- 
based approaches. Usually these are non-self-consistent 
techniques that approach the problem from the point of 
view of transmission and reflection of Bloch waves moving 
through the device. We found that because the structure 
in the leads is frozen in beginning with A'^ = 1 and be- 
cause the exponential decay lengths are also determined 
from iV = 1, if one knew those results and plugged them 
into the Landauer approach, one should be able to cal- 
culate accurate properties; i.e. the self consistency is 
needed for each nanostructure, but the self-consistency 
hardly changes with the thickness of the barrier. Hence 
a phenomenological approach that adjusts the properties 
of the barrier height to produce the required behavior, 
may work well, even for strongly correlated systems; of 
course, the many-body theory is the only way to deter- 
mine the precise structure needed via its self-consistent 
solution (i.e. it requires no fitting). 

There are a number of important effects that we have 



not discussed here, which play roles in the transport 
through nanostructures. We did not attempt to include 
them in this first, simplest problem that we tackled. The 
first one is the issue of charge reorganization around the 
interface. If the chemical potentials of the leads and the 
barriers are different, electrons will spill from one plane 
to the another until a screened dipole layer is formed, and 
a constant electrochemical potential is found throughout 
the deviceSS. Such effects can have dramatic results if one 
or more of the materials is a correlated insulator, since 
the inhomogeneous doping of the system can transform 
part of it from insulating to metallic. This is believed 
to occur in grain boundaries in high temperature su- 
perconducting tapes and wires2i, and in insulator-based 
nanostructures^^'^'^. Second, calculations should be per- 
formed off of half filling, where the thermal evolution 
of the chemical potential, will likely undergo some tem- 
perature dependence so the charge rearrangement can 
vary with temperature in the system. Third, we should 
calculate the thermal transport effects. Since these cal- 
culations require particle-hole asymmetry, we will have 
the chemical potential evolution and the charge reorga- 
nizations to deal with as well. Fourth, one can include 
ordered phase effects at the mean-field level easily, as in 
a superconductor for a Josephson junctionSi, or in a fer- 
romagnet for a spintronics device. Fifth, it will be useful 
to determine the capacitance of a nanostructure, since 
the capacitance is often important in determining the 
switching speed of a device; it can be calculated with a 
linear-response formalism as well. Finally, we also should 
look into nonequilibrium effects, especially the nonlinear 
response of a current-voltage curve. It is our plan to 
investigate these complications in the future. 
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